function [betas,gamma,c,fX,eVector]=estVlstar(y,x,q,c)

N=size(x,2);

% Compute starting values for gamma and c using grid search
[betas,gamma,c] = startvalVlstar(y,x,q,c);

options = optimset('MaxFunEvals',1000,'TolFun',1e-30, 'TolX',1e-30,'Display','off');

gamma = lsqnonlin(@(gamma)eVlstar(gamma,y,x,q,c),gamma,[],[],options);
 
%computing betas using OLS given gamma and c          
fX  = 1./(1+exp(-exp(gamma)/std(q)*(q-c)));

rec = 1-fX;

T=size(y,1);
N=size(y,2);

betas=[];
e=[];
for i=1:N;

    y_i=y(:,i);

    x_i=[ones(T,1) x(:,i)];

    x_i=[x_i.*(rec*ones(1,2)) x_i.*(1-(rec(:,1)*ones(1,2)))];

    betas_i=x_i\y_i;

    betas=[ betas; betas_i ];

    e_i=y_i-x_i*betas_i;
    e=[ e e_i ];

end

eVector=reshape(e,[T*N,1]); 


